/*
*    OpenCL program to compute convex hull of 2D points using gift wrapping method.
*	 It is compiled on Windows system but It can be easily ported to Linux system.
*
*    Copyright (C) 2013, Ahmet Artu Yildirim
*
*    This program is free software: you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation, either version 3 of the License, or
*    (at your option) any later version.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*    GNU General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <windows.h>
#include <CL/opencl.h>

typedef struct _point {
	float x;
	float y;
} point;

LARGE_INTEGER frequency;        // ticks per second
    LARGE_INTEGER t1, t2;           // ticks
    double elapsedTime;

const point POINTNULL = {FLT_MIN, FLT_MIN};
const int WIDTH = 20000;
const int HEIGHT = 20000;

point* quickhull(point* s, int s_length);
point* giftwrapping(point* s, int s_length);
void createsvg(point* s, int s_length, point* hull);

point* quickhull(point* S, int lengthofS) {
	int i, j, res;
	int s_length = lengthofS;
	int anslength = 0;
	int ischanged = 1;
	int* label = (int*) malloc(s_length * sizeof (int));
	int* state = (int*) malloc(s_length * sizeof (int));
	float* dist = (float*) malloc(s_length * sizeof (float));
	float* maxdist = (float*) malloc(s_length * sizeof (float));
	point* maxdistpoints = (point*) malloc(s_length * sizeof (point));
	point* ans = (point*) malloc(s_length * sizeof (point));
	point* P = (point*) malloc(s_length * sizeof (point));
	int mdlength = 0;

	point min = {FLT_MAX, FLT_MAX};
	point max = {FLT_MIN, FLT_MIN};

	/* initialize label, ans and find min max points*/
	for (i = 0; i < s_length; i++) {
		label[i] = 0;
		P[i] = S[i];
		ans[i] = POINTNULL;

		if (S[i].x < min.x)
			min = S[i];

		if (S[i].x > max.x)
			max = S[i];
	}

	ans[anslength++] = min;
	ans[anslength++] = max;

	printf("MIN: X: %lg Y: %lg\n", min.x, min.y);
	printf("MAX: X: %lg Y: %lg\n", max.x, max.y);

	while (ischanged) {
		for (i = 0; i < s_length; i++) {
			point cur = P[i];
			int curlabel = label[i];
			point l = ans[curlabel];
			point r = ans[curlabel + 1];
			float s = (r.y - l.y) / (r.x - l.x);
			float d = s * (cur.x - l.x) - cur.y + r.y;
			dist[i] = d;

			printf("X: %lg Y: %lg Dist: %lg Label: %d\n", cur.x, cur.y, d, label[i]);
		}

		/* clear max dist and state */
		for (i = 0; i < anslength; i++) {
			maxdist[i] = 0.0f;
			state[i] = 0;
		}

		ischanged = 0;

		for (i = 0; i < s_length; i++) {
			int curlabel = label[i];
			point cur = P[i];

			if (maxdist[curlabel] > dist[i]) {
				maxdistpoints[curlabel] = cur;
				maxdist[curlabel] = dist[i];
				state[curlabel] = 1;
				ischanged = 1;

				printf("maxdist %d: X: %lg Y: %lg\n", j, P[i].x, P[i].y);
			}
		}

		res = 0;

		for (i = 0; i < anslength; i++) {
			res = res + state[i];
			state[i] = res - state[i];
		}

		/* copy points which has distance greater than zero*/
		j = 0;
		for (i = 0; i < s_length; i++) {
			if (dist[i] >= 0.0f) {
				P[j++] = P[i];
			}
		}
		s_length = j;

		for (i = 0; i < s_length; i++) {
			int curlabel = label[i] + state[label[i]];

			if (P[i].x >= ans[curlabel + 1].x) {
				label[i] = curlabel + 1;
			} else {
				label[i] = curlabel;
			}

			//printf ("%d Point: %d\n", i, label[i]);
		}

		j = 0;
		//printf("Max dist length: %d\n", mdlength);

		for (i = 0; i < anslength; i++) {
			ans[j] = maxdistpoints[j];

			//printf ("Result %d: X: %lg Y: %lg\n", j, ans[i].x, ans[i].y);
		}

		// printf ("Result %d: X: %lg Y: %lg\n", j, ans[i].x, ans[i].y);

		if (anslength != lengthofS)
			ans[anslength - 1] = POINTNULL;
	}

	for (i = 0; i < s_length; i++) {
		point cur = P[i];
		printf("X: %lg Y: %lg label: %d\n", cur.x, cur.y, label[i]);
	}

	return ans;
}

float* giftwrapping(float* s, int s_length) {
	if (s_length < 3) {
		printf("There must be 3 points at least!");
		exit(-1);
	}

	int i;
	float* ans = (float*) malloc((s_length + 2) * 2 * sizeof (float));
	int ansindex = 0;

	/* find left most point and right most point*/
	int min = 0;
	int max = 0;

	for (i = 1; i < s_length; i++) {
		if (s[i * 2] == s[min * 2]) {
			if (s[i * 2 + 1] < s[min * 2 + 1])
				min = i;
		} else if (s[i * 2] < s[min * 2])
			min = i;

		if (s[i * 2] == s[max * 2]) {
			if (s[i * 2 + 1] < s[max * 2 + 1])
				max = i;
		} else if (s[i * 2] > s[max * 2])
			max = i;
	}

	int smallest;
	int current = min;
	float line1Start[2];
	float line1End[2];
	float line2Start[2];
	float line2End[2];
	float smallestangle;
	int reverse = 0;

	QueryPerformanceCounter(&t1);

	do {
		ans[ansindex * 2] = s[current * 2];
		ans[ansindex * 2 + 1] = s[current * 2 + 1];
		ansindex++;

		smallestangle = -360.0f;

		for (i = 0; i < s_length; i++) {
			if (current == i)
				continue;

			if (reverse) {
				line1Start[0] = WIDTH - s[current * 2];
				line1Start[1] = HEIGHT - s[current * 2 + 1];
				line1End[0] = line1Start[0] + WIDTH;
				line1End[1] = line1Start[1];
				line2Start[0] = line1Start[0];
				line2Start[1] = line1Start[1];
				line2End[0] = WIDTH - s[i * 2];
				line2End[1] = HEIGHT - s[i * 2 + 1];
			} else {
				line1Start[0] = s[current * 2];
				line1Start[1] = s[current * 2 + 1];
				line1End[0] = line1Start[0] + WIDTH;
				line1End[1] = line1Start[1];
				line2Start[0] = line1Start[0];
				line2Start[1] = line1Start[1];
				line2End[0] = s[i * 2];
				line2End[1] = s[i * 2 + 1];
			}

			float atanA = atan2(line1End[0] - line1Start[0], line1End[1] - line1Start[1]);
			float atanB = atan2(line2End[0] - line2Start[0], line2End[1] - line2Start[1]);

			float cangle = (atanB - atanA) * 57.2957795786f;

			if (cangle > smallestangle) {
				smallestangle = cangle;
				smallest = i;
			}
		}
		current = smallest;

		if (current == max)
			reverse = 1;

	} while (current != min);

	QueryPerformanceCounter(&t2);

	ans[ansindex * 2] = ans[0];
	ans[ansindex * 2 + 1] = ans[1];
	ans[ansindex * 2 + 2] = FLT_MIN;

	return ans;
}

void createsvg(float* s, int s_length, float* hull) {
	FILE *ofp;

	ofp = fopen("out.svg", "w");

	if (ofp == NULL) {
		fprintf(stderr, "Can't open output file!");
		exit(1);
	}

	int i;
	const int sfactor = 10;

	/* write svg header*/
	fprintf(ofp, "<?xml version=\"1.0\" standalone=\"no\" ?>\n");
	fprintf(ofp, "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.0//EN\"\n");
	fprintf(ofp, "\t\"http://www.w3.org/TR/2001/REC-SVG-20010904/DTD/svg10.dtd\">\n\n");

	/* fill svg*/
	fprintf(ofp, "<svg width=\"%d\" height=\"%d\" xmlns=\"http://www.w3.org/2000/svg\">\n", WIDTH * sfactor, HEIGHT * sfactor);

	/* draw points as circle */
	for (i = 0; i < s_length; i++) {
		fprintf(ofp, "<circle cx=\"%d\" cy=\"%d\" r=\"%d\" stroke=\"red\" stroke-width=\"1\" fill=\"red\"/>\n", (int) (s[i * 2] * sfactor), (int) (s[i * 2 + 1] * sfactor), sfactor / 2);
	}

	/* draw lines */
	i = 0;
	while (1) {
		float p1x = hull[i * 2];
		if (p1x == FLT_MIN)
			break;

		float p1y = hull[i * 2 + 1];
		float p2x = hull[i * 2 + 2];

		if (p2x == FLT_MIN)
			break;

		float p2y = hull[i * 2 + 3];

		fprintf(ofp, "<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" style=\"stroke:rgb(255,0,0);stroke-width:2\"/>\n",
			(int) (p1x * sfactor), (int) (p1y * sfactor), (int) (p2x * sfactor), (int) (p2y * sfactor));
		i++;
	}

	fprintf(ofp, "</svg>");
	fclose(ofp);
}

const char *kernelSource =													   "\n" \
	"__kernel void findMinAngle(__global float *s,									\n" \
	"	__global int *current,													    \n" \
	"	__global int *reverse,														\n" \
	"	const unsigned int n,														\n" \
	"	const unsigned int width,													\n" \
	"	const unsigned int height,													\n" \
	"	__global float *out_smallest)												\n" \
	" {																				\n" \
	"	//Get our global thread ID													\n" \
	"	int id = get_global_id(0);int local_index = get_local_id(0);													\n" \
	"	__local float localangles[512];	float cangle = -360.0f;											\n" \
	"	//Make sure we do not go out of bounds										\n" \
	"	if (id < n && id != *current)												\n" \
	"	{																			\n" \
	"		float line1Startx;														\n" \
	"		float line1Starty;														\n" \
	"		float line1Endx;														\n" \
	"		float line1Endy;														\n" \
	"		float line2Startx;														\n" \
	"		float line2Starty;														\n" \
	"		float line2Endx;														\n" \
	"		float line2Endy;														\n" \
	"																				\n" \
	"		if (*reverse) {															\n" \
	"			line1Startx = width - s[*current * 2];								\n" \
	"			line1Starty = height - s[*current * 2 + 1];							\n" \
	"			line1Endx = line1Startx + width;									\n" \
	"			line1Endy = line1Starty;											\n" \
	"			line2Startx = line1Startx;											\n" \
	"			line2Starty = line1Starty;											\n" \
	"			line2Endx = width - s[id * 2];										\n" \
	"			line2Endy = height - s[id * 2 + 1];									\n" \
	"		} else {																\n" \
	"			line1Startx = s[*current * 2];										\n" \
	"			line1Starty = s[*current * 2 + 1];									\n" \
	"			line1Endx = line1Startx + width;									\n" \
	"			line1Endy = line1Starty;											\n" \
	"			line2Startx = line1Startx;											\n" \
	"			line2Starty = line1Starty;											\n" \
	"			line2Endx = s[id * 2];												\n" \
	"			line2Endy = s[id * 2 + 1];											\n" \
	"		}																		\n" \
	"																				\n" \
	"		float atanA = atan2(line1Endx - line1Startx, line1Endy - line1Starty);	\n" \
	"		float atanB = atan2(line2Endx - line2Startx, line2Endy - line2Starty);	\n" \
	"																				\n" \
	"		cangle = (atanB - atanA) * 57.2957795786f;						\n" \
	"																				\n" \
	"		// Perform parallel reduction											\n" \
	"												\n" \
	"  }																		\n" \
	"		localangles[local_index * 2] = cangle;									\n" \
	"		localangles[local_index * 2 + 1] = id;									\n" \
	"		barrier(CLK_LOCAL_MEM_FENCE);											\n" \
	"		for(int offset = get_local_size(0) / 2; offset > 0; offset = offset / 2)\n" \
	"		{																		\n" \
	"			if (local_index < offset)											\n" \
	"			{																	\n" \
	"	if (id < n && id != *current)												\n" \
	"	{																			\n" \
	"				float other = localangles[(local_index + offset) * 2];			\n" \
	"				float mine = localangles[local_index * 2];						\n" \
	"				if (mine < other)												\n" \
	"				{																\n" \
	"					localangles[local_index * 2] = other;						\n" \
	"					localangles[local_index * 2 + 1] = localangles[(local_index + offset) * 2 + 1];				\n" \
	"				}																\n" \
	"  }																		\n" \
	"			}																	\n" \
	"			barrier(CLK_LOCAL_MEM_FENCE);										\n" \
	"		}																		\n" \
	"		if (local_index == 0)													\n" \
	"		{																		\n" \
	"			out_smallest[get_group_id(0) * 2] = localangles[0];					\n" \
	"			out_smallest[get_group_id(0) * 2 + 1] = localangles[1];				\n" \
	"		}																		\n" \
	"																				\n" \
	"}																				\n" \
	"\n" ;


float* giftwrapping_opencl(float* s, int s_length) {
	if (s_length < 3) {
		printf("There must be 3 points at least!");
		exit(-1);
	}

	int i;
	float* ans = (float*) malloc((s_length + 2) * 2 * sizeof (float));
	int ansindex = 0;

	/* find left most point and right most point*/
	int min = 0;
	int max = 0;

	for (i = 1; i < s_length; i++) {
		if (s[i * 2] == s[min * 2]) {
			if (s[i * 2 + 1] < s[min * 2 + 1])
				min = i;
		} else if (s[i * 2] < s[min * 2])
			min = i;

		if (s[i * 2] == s[max * 2]) {
			if (s[i * 2 + 1] < s[max * 2 + 1])
				max = i;
		} else if (s[i * 2] > s[max * 2])
			max = i;
	}

	int smallest;
	int current = min;
	float line1Start[2];
	float line1End[2];
	float line2Start[2];
	float line2End[2];
	float smallestangle;
	int reverse = 0;

	// Device input buffers
	cl_mem d_s;
	cl_mem d_current;
	cl_mem d_reverse;
	// Device output buffer
	cl_mem d_smallestresult;

	cl_platform_id cpPlatform;        // OpenCL platform
	cl_device_id device_id;           // device ID
	cl_context context;               // context
	cl_command_queue queue;           // command queue
	cl_program program;               // program
	cl_kernel kernel;                 // kernel
	cl_device_id device_ids[2];

	size_t globalSize, localSize;
	int numberofworkgroups;
	cl_int err;

	// Bind to platform
	cl_uint numberofplatforms = 0;
	err = clGetPlatformIDs(1, NULL, &numberofplatforms);

	// assumed there are two opencl device. on my machine first one is GPU and the second one is CPU.
	cl_platform_id platformids[2];
	bool isgpu = true;
	if (isgpu) {
		err = clGetPlatformIDs(2, platformids, NULL);

		cpPlatform = platformids[0];
		if (CL_SUCCESS != err) {
			printf("clGetPlatformIDs failed\n");
			exit(-1);
		}
		// Get ID for the device
		err = clGetDeviceIDs(cpPlatform, CL_DEVICE_TYPE_GPU, 1, &device_id, NULL);
	} else {
			err = clGetPlatformIDs(2, platformids, NULL);

		cpPlatform = platformids[1];

			if (CL_SUCCESS != err) {
				printf("clGetPlatformIDs failed\n");
				exit(-1);
			}
			// Get ID for the device
			err = clGetDeviceIDs(cpPlatform, CL_DEVICE_TYPE_CPU, 1, &device_id, NULL);
	}

	if (CL_SUCCESS != err) {
		if (err == CL_INVALID_PLATFORM)
		{
			printf("clGetDeviceIDs failed: CL_INVALID_PLATFORM\n");
		} else if (err == CL_INVALID_DEVICE_TYPE )
		{
			printf("clGetDeviceIDs failed: CL_INVALID_DEVICE_TYPE \n");
		} else if (err == CL_INVALID_VALUE  )
		{
			printf("clGetDeviceIDs failed: CL_INVALID_VALUE  \n");
		}else if (err == CL_DEVICE_NOT_FOUND   )
		{
			printf("clGetDeviceIDs failed: CL_DEVICE_NOT_FOUND   \n");
		}else if (err == CL_OUT_OF_RESOURCES    )
		{
			printf("clGetDeviceIDs failed: CL_OUT_OF_RESOURCES    \n");
		}else if (err == CL_OUT_OF_HOST_MEMORY     )
		{
			printf("clGetDeviceIDs failed: CL_OUT_OF_HOST_MEMORY     \n");
		} else {
			printf("clGetDeviceIDs failed\n");
		}
		exit(-1);
	}

	// Get max dimensions supported by the device
	cl_uint maxdimensionssupported;
	clGetDeviceInfo(device_id, CL_DEVICE_MAX_WORK_ITEM_DIMENSIONS, sizeof(cl_uint), &maxdimensionssupported, NULL);

	// Get max work item sizes for all dimensions
	size_t* maxworkitemsizes = (size_t*) malloc(maxdimensionssupported * sizeof(size_t));
	clGetDeviceInfo(device_id, CL_DEVICE_MAX_WORK_ITEM_SIZES, maxdimensionssupported * sizeof(size_t), maxworkitemsizes, NULL);

	localSize = 256;//maxworkitemsizes[0];

	free (maxworkitemsizes);

	// Number of total work items - localSize must be devisor
	globalSize = ceil(s_length/(float)localSize)*localSize;

	numberofworkgroups = globalSize / localSize;

	// Create a context
	context = clCreateContext(0, 1, &device_id, NULL, NULL, &err);

	// Create a command queue
	queue = clCreateCommandQueue(context, device_id, 0, &err);

	// Create the compute program from the source buffer
	program = clCreateProgramWithSource(context, 1,
		(const char **) & kernelSource, NULL, &err);

	// Build the program executable
	cl_int ret = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);

	if(ret != CL_SUCCESS)
	{
		printf("clBuildProgram failed: %d\n", ret);
		char buf[0x10000];
		clGetProgramBuildInfo( program,
			device_id,
			CL_PROGRAM_BUILD_LOG,
			0x10000,
			buf,
			NULL);
		printf("\n%s\n", buf);
		exit(-1);
	}

	// Create the compute kernel in the program we wish to run
	kernel = clCreateKernel(program, "findMinAngle", &err);

	QueryPerformanceCounter(&t1);

	// Create the input and output arrays in device memory for our calculation
	d_s = clCreateBuffer(context, CL_MEM_READ_ONLY, s_length * 2 * sizeof(float), NULL, NULL);
	d_current = clCreateBuffer(context, CL_MEM_READ_ONLY, sizeof(int), NULL, NULL);
	d_reverse = clCreateBuffer(context, CL_MEM_READ_ONLY, sizeof(int), NULL, NULL);

	d_smallestresult = clCreateBuffer(context, CL_MEM_WRITE_ONLY, numberofworkgroups * 2 *sizeof(float), NULL, NULL);

	float * smallestresult = (float*) malloc (numberofworkgroups * 2 *sizeof(float));

	// Set the arguments to our compute kernel
	err  = clSetKernelArg(kernel, 0, sizeof(cl_mem), &d_s);
	err |= clSetKernelArg(kernel, 1, sizeof(cl_mem), &d_current);
	err |= clSetKernelArg(kernel, 2, sizeof(cl_mem), &d_reverse);
	err |= clSetKernelArg(kernel, 3, sizeof(unsigned int), &s_length);
	err |= clSetKernelArg(kernel, 4, sizeof(unsigned int), &WIDTH);
	err |= clSetKernelArg(kernel, 5, sizeof(unsigned int), &HEIGHT);
	err |= clSetKernelArg(kernel, 6, sizeof(cl_mem), &d_smallestresult);

	// Write input buffer, current id and reverse status in device memory
	err = clEnqueueWriteBuffer(queue, d_s, CL_TRUE, 0,
		s_length * 2 * sizeof(float), s, 0, NULL, NULL);
	err |= clEnqueueWriteBuffer(queue, d_reverse, CL_TRUE, 0,
		sizeof(int), &reverse, 0, NULL, NULL);

	do {
		ans[ansindex * 2] = s[current * 2];
		ans[ansindex * 2 + 1] = s[current * 2 + 1];
		ansindex++;

		err |= clEnqueueWriteBuffer(queue, d_current, CL_TRUE, 0,
			sizeof(int), &current, 0, NULL, NULL);

		if (current == max)
		{
			reverse = 1;
			err |= clEnqueueWriteBuffer(queue, d_reverse, CL_TRUE, 0,
				sizeof(int), &reverse, 0, NULL, NULL);
		}

		// Execute the kernel over the entire range of the data set
		err = clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &globalSize, &localSize,
			0, NULL, NULL);

		// end of kernel invocation
		// Wait for the command queue to get serviced before reading back results
		clFinish(queue);

		// Read the element id with smallest angle from the device
		clEnqueueReadBuffer(queue, d_smallestresult, CL_TRUE, 0,
			numberofworkgroups * 2 * sizeof(float), smallestresult, 0, NULL, NULL );

		// Reduce group results
		smallestangle = smallestresult[0];
		smallest = (int) smallestresult[1];

		for (i = 1; i < numberofworkgroups; i++)
		{
			if (smallestresult[i * 2] > smallestangle)
			{
				smallestangle = smallestresult[i * 2];
				smallest = (int) smallestresult[i * 2 + 1];
			}
		}

		current = smallest;

	} while (current != min);

	QueryPerformanceCounter(&t2);

	ans[ansindex * 2] = ans[0];
	ans[ansindex * 2 + 1] = ans[1];
	ans[ansindex * 2 + 2] = FLT_MIN;

	// release OpenCL resources
	clReleaseMemObject(d_s);
	clReleaseMemObject(d_current);
	clReleaseMemObject(d_reverse);
	clReleaseMemObject(d_smallestresult);
	clReleaseProgram(program);
	clReleaseKernel(kernel);
	clReleaseCommandQueue(queue);
	clReleaseContext(context);

	// release host memory
	free (smallestresult);

	return ans;
}

/*
* 
*/
int main(int argc, char** argv) {
	QueryPerformanceFrequency(&frequency);

	int lengthofs = 262144;
	int i;

	float* s = (float*) malloc(lengthofs * 2 * sizeof (float));
	/* generate random points */
	for (i = 0; i < lengthofs; i++) {
		s[i * 2] = (float) (rand() % WIDTH);
		s[i * 2 + 1] = (float) (rand() % HEIGHT);
	}

	/*
	s[0].x = 5;
	s[0].y = 15;
	s[1].x = 15;
	s[1].y = 5;
	s[2].x = 20;
	s[2].y = 10;
	s[3].x = 25;
	s[3].y = 5;
	s[4].x = 30;
	s[4].y = 15;
	s[5].x = 15;
	s[5].y = 20;
	s[6].x = 10;
	s[6].y = 10;
	s[7].x = 27;
	s[7].y = 7;*/

	float* ans = giftwrapping(s, lengthofs); //giftwrapping_opencl(s, lengthofs);

	if (ans != NULL);
		createsvg(s, lengthofs, ans);

	elapsedTime = (t2.QuadPart - t1.QuadPart) * 1000.0 / frequency.QuadPart;
	printf("elapsedTime: %f\n", elapsedTime);

	return (EXIT_SUCCESS);
}

